!>\brief
!!
!! Poisson equation.
!!
!! \n
!!
!! In this module we provide some Poisson solvers. The notation and
!! the main definitions are those introduced in \c mod_vp_ode.
!!
!! \section DGclassic The classical DG approach
!!
!! We consider here the classical DG formulation, where both primal
!! and dual variables belong to the same polynomial space.
!!
!! In matrix form, the Poisson problem for the computation of the
!! electric field \f$\underline{E}\f$ reads
!! \f{displaymath}{
!!  \left[\begin{array}{cc}
!!   M & D^{LDG,T} \\
!!   -D^{LDG} & C
!!  \end{array}\right]
!!  \left[\begin{array}{c}
!!   \underline{E} \\
!!   \underline{\Phi}
!!  \end{array}\right] =
!!  \left[\begin{array}{c}
!!   0 \\
!!   \underline{f}
!!  \end{array}\right].
!! \f}
!! The computation of the mass matrix is trivial (setting
!! \f$\underline{q}_{k\underline{i}} =
!! \phi_{\underline{i}}\underline{e}_{k}\f$):
!! \f{displaymath}{
!!   M_{k\underline{i}\,h\underline{j}} =
!!   \int_{K_x}\underline{q}_{k\underline{i}}\cdot
!!   \underline{q}_{h\underline{j}}\,dx = \sum_{\underline{l}}
!!   w_{\underline{l}} J_x
!!   \phi_{\underline{i}}(\underline{x}_{\underline{l}}^G)
!!   \phi_{\underline{j}}(\underline{x}_{\underline{l}}^G)\,
!!   \delta_{kh} =
!!   J_xw_{\underline{i}}\delta_{\underline{i}\underline{j}} \,
!!   \delta_{kh}.
!! \f}
!! Concerning the differentiation matrix:
!! \f{displaymath}{
!!   D^{LDG}_{\underline{i}\,h\underline{j}} = - \int_{K_x}
!!  \underline{q}_{h\underline{j}}\cdot\nabla\phi_{\underline{i}} \, dx
!!  +\int_{\partial K_x}\left( \{\underline{q}_{h\underline{j}}\} -
!!  \underline{c}_{12}\llbracket\underline{q}_{h\underline{j}}\rrbracket
!!  \right)\cdot\underline{n}\phi_{\underline{i}}\,d\sigma
!! \f}
!! The firs term can be computed as a special case of the volume term
!! in the Vlasov equation; the result is zero unless
!! \f$\underline{i}_{\backslash h} =
!! \underline{j}_{\backslash h}\f$, and the nonzero terms are
!! \f{displaymath}{
!!   - \int_{K_x}
!!   \underline{q}_{h\underline{j}}\cdot\nabla\phi_{\underline{i}} \,
!!   dx = - \frac{J_xw_{\underline{i}}}{J_hw_{\underline{i}(h)}}
!!   D_{\underline{i}(h)\underline{j}(h)}.
!! \f}
!! Concerning the boundary term, first of all we notice that this is
!! nonzero only if \f$\underline{i}\in s_x\f$ for some side \f$s_x\f$.
!! Then, in such case we have
!! \f{displaymath}{
!!  \int_{s_x}\left( \{\underline{q}_{h\underline{j}}\} -
!!  \underline{c}_{12}\llbracket\underline{q}_{h\underline{j}}\rrbracket
!!  \right)\cdot\underline{n}_{\underline{i}}
!!  \phi_{\underline{i}}\,d\sigma
!!  = w_{\underline{i}_{\backslash n^s}}|s_x|\left( \frac{1}{2} -
!!  \underline{c}_{12}\cdot\underline{n}_{\underline{j}} \right)
!!  \phi_{\underline{j}}(\underline{x}^G_{\underline{i}})\,
!!  (\underline{n}_{\underline{i}})_h,
!! \f}
!! where a subscript is added to the side normal to indicate that
!! \f$\underline{n}_{\underline{i}}\f$ is directed outward with
!! respect to the element containing
!! \f$\underline{x}^G_{\underline{i}}\f$. Since for each quadrature
!! node on \f$s_x\f$ there are two possible choices for
!! \f$\underline{i}\f$ and \f$\underline{j}\f$, a two-by-two block
!! results.
!! Finally, one has to consider the stabilization matrix
!! \f{displaymath}{
!!   C_{\underline{i}\underline{j}} = \int_{\partial K_x} c_{11}
!!   \llbracket \phi_{\underline{j}} \rrbracket \cdot\underline{n}
!!   \phi_{\underline{j}}\,d\sigma
!!  = w_{\underline{i}_{\backslash n^s}}|s_x|
!!  \left(
!!  \underline{n}_{\underline{j}}\cdot\underline{n}_{\underline{i}}
!!  \right) c_{11}
!!  \phi_{\underline{j}}(\underline{x}^G_{\underline{i}}).
!! \f}
!!
!! \subsection impdetails Implementation details
!!
!! To ensure that the average of \f$\Phi\f$ is yero, we introduce a
!! Lagrange multiplier. This means that the system is redefined as
!! \f{displaymath}{
!!  \left[\begin{array}{ccc}
!!   M & D^{LDG,T} \\
!!   -D^{LDG} & C & L^T \\
!!       & L &
!!  \end{array}\right]
!!  \left[\begin{array}{c}
!!   \underline{E} \\
!!   \underline{\Phi} \\
!!   \lambda
!!  \end{array}\right] =
!!  \left[\begin{array}{c}
!!   0 \\
!!   \underline{f}\\
!!   0
!!  \end{array}\right]
!! \f}
!! with
!! \f{displaymath}{
!!  L_{\underline{j}} = \int_{K_x}\phi_{\underline{j}}\,dx =
!!  J_xw_{\underline{j}}.
!! \f}
!! This can be recast in the original form by setting
!! \f{displaymath}{
!!  \tilde{D}^{LDG} = \left[ \begin{array}{c}
!!   D^{LDG} \\ 0
!!  \end{array}\right], \qquad
!!  \tilde{C} = \left[ \begin{array}{cc}
!!   C & L^T \\ L & 0
!!  \end{array}\right], \qquad
!!  \tilde{\underline{\Phi}} = \left[ \begin{array}{c}
!!   \underline{\Phi} \\ \lambda
!!  \end{array}\right], \qquad
!!  \tilde{\underline{f}} = \left[ \begin{array}{c}
!!   \underline{f} \\ 0
!!  \end{array}\right],
!! \f}
!! so that the reduced system reads
!! \f{displaymath}{
!!  \left( \tilde{D}^{LDG}\,M^{-1}\,\tilde{D}^{LDG,T} +
!!  \tilde{C}\right) \tilde{\underline{\Phi}} = \tilde{\underline{f}}.
!! \f}
!! The electric field can then be computed as
!! \f{displaymath}{
!!  \underline{E} = -M^{-1}\tilde{D}^{LDG,T}\tilde{\underline{\Phi}}.
!! \f}
!!
!! \section DGwRT DG method with Raviart-Thomas basis functions
!!
!! We consider here an approach obtained modifying the \ref DGclassic
!! "classical DG approach" in the following way:
!! <ul>
!!  <li> the flux unknown belongs to the Raviart-Thomas space
!!  <li> suitable quadrature formulas are used for each term of the
!!  equations.
!! </ul>
!! More precisely, we consider the problem
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!  \displaystyle
!!  \int_{\underline{Q}_1,K_x}\underline{E}\cdot\underline{z}\,dx
!!  -\int_{\underline{Q}_1,K_x}\nabla\Phi\cdot\underline{z}\,dx
!!  -\int_{Q_s,\partial
!!  K_x}\left(\widehat{\Phi}-\Phi\right)\underline{n}\cdot\underline{z}
!!  \, d\sigma & = & 0 \\[4mm] \displaystyle
!!  \int_{\underline{Q}_1,K_x}\underline{E}\cdot\nabla p\,dx -
!!  \int_{Q_s,\partial
!!  K_x}\underline{n}\cdot\widehat{\underline{E}}\,p\,d\sigma & = &
!!  \displaystyle -\int_{Q_3,K_x} (1-\rho)\,p\,dx,
!!  \end{array}
!! \f}
!! where \f$Q_3\f$ and \f$Q_s\f$ are the standard quadrature formulas
!! used in the spectral DG method, while \f$\underline{Q}_1\f$ is a
!! quadrature formula of increased accuracy required to integrate the
!! Raviart-Thomas polynomials. Starting with the mass matrix, we have
!! \f{displaymath}{
!!   M_{k\underline{i}^{k+}\,h\underline{j}^{h+}} =
!!  \int_{\underline{Q}_1,K_x}\underline{q}_{k\underline{i}^{k+}}\cdot \,
!!   \underline{q}_{h\underline{j}^{h+}}\,dx =
!!   \delta_{kh}
!!  \int_{\underline{Q}_1,K_x}\phi_{\underline{i}^{k+}}
!!   \phi_{\underline{j}^{k+}}\,dx =
!!   \delta_{kh}
!!   \sum_{\underline{l}^{k+}}
!!   w_{\underline{l}^{k+}} J_x
!!   \phi_{\underline{i}^{k+}}(\underline{x}_{\underline{l}^{k+}}^G)
!!   \phi_{\underline{j}^{k+}}(\underline{x}_{\underline{l}^{k+}}^G).
!! \f}
!! In this last expression, we have introduced a new multiindex
!! \f$\underline{l}^{k+}\f$, such that its components are
!! \f{displaymath}{
!!  \underline{l}^{k+}(i) = \left\{
!!   \begin{array}{cl}
!!    1,\ldots,r+1, & i\neq k \\
!!    1,\ldots,r+2, & i= k
!!   \end{array} \right.
!! \f}
!! having denoted by \f$r\f$ the polynomial degree. This extended
!! multiindex enumerates the degrees of freedom of the Raviart-Thomas
!! space of order \f$r\f$. The multidimensional Gaussian nodes and
!! weights are defined as in the standard DG case, except that the
!! one-dimensional nodes and weights are taken from <em>two</em>
!! quadrature formulas with \f$r+1\f$ and \f$r+2\f$ nodes,
!! respectively. The same consideration holds for the associated basis
!! functions, which are build as tensor products of the
!! one-dimensional basis functions of <em>two</em> basis. The final
!! result is
!! \f{displaymath}{
!!   M_{k\underline{i}^{k+}\,h\underline{j}^{h+}} =
!!   J_xw_{\underline{i}^{k+}}
!!   \delta_{\underline{i}^{k+}\underline{j}^{k+}} \, \delta_{kh},
!! \f}
!! which is formally identical to the classical DG case.
!! \note This implies that, when integrating a scalar product
!! \f$\underline{a}\cdot\underline{b} = \sum_{i=1}^d a_ib_i\f$, each
!! term \f$a_ib_i\f$ is integrated with a specific quadrature formula,
!! using nodes and weights with index \f$\underline{l}^{i+}\f$. In
!! fact, this is the definition of the quadrature formula
!! \f$\underline{Q}_1\f$, and clarifies the use of the underlined
!! notation, in order to emphasize that we have a vector of quadrature
!! formulas.
!!
!! \note A critical point in obtaining a diagonal mass matrix is that
!! basis functions with index \f$\underline{i}^{k+}\f$ are never
!! multiplied with basis functions with index
!! \f$\underline{j}^{h+}\f$ with \f$h\neq k\f$.
!!
!! The situation is not as simple for the differentiation matrix,
!! because now one has to multiply basis functions belonging to
!! different basis. To see this, let us consider
!! \f{displaymath}{
!!   D^{LDG}_{\underline{i}\,h\underline{j}^{h+}} = -
!!   \int_{\underline{Q}_1,K_x}
!!   \underline{q}_{h\underline{j}^{h+}}\cdot\nabla\phi_{\underline{i}}
!!   \, dx +\int_{Q_s,\partial K_x}\left(
!!   \{\underline{q}_{h\underline{j}^{h+}}\} - \underline{c}_{12}
!!   \llbracket\underline{q}_{h\underline{j}^{h+}}\rrbracket
!!   \right)\cdot\underline{n}\phi_{\underline{i}}\,d\sigma.
!! \f}
!! Taking the volume integral we have
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!  \displaystyle
!!    - \int_{\underline{Q}_1,K_x}
!!    \underline{q}_{h\underline{j}^{h+}}\cdot\nabla\phi_{\underline{i}}
!!    \, dx & = &
!!  \displaystyle
!!  - \sum_{k=1}^d
!!    \int_{\underline{Q}_1(k),K_x}
!!    \left(\underline{q}_{h\underline{j}^{h+}}\right)_k
!!    \frac{\partial\phi_{\underline{i}}}{\partial x_k}
!!    \, dx \\[5mm]
!! & = &
!!  \displaystyle
!!  - \int_{\underline{Q}_1(h),K_x}
!!    \phi_{\underline{j}^{h+}}
!!    \frac{\partial\phi_{\underline{i}}}{\partial x_h}
!!    \, dx \\[5mm]
!! & = &
!!  \displaystyle
!!  - J_x \sum_{\underline{l}^{h+}} w_{\underline{l}^{h+}}\,
!!    \phi_{\underline{j}^{h+}}(\underline{x}^G_{\underline{l}^{h+}})\,
!!    \frac{\partial\phi_{\underline{i}}}{\partial x_h}
!!    (\underline{x}^G_{\underline{l}^{h+}}) \\[5mm]
!! & = &
!!  \displaystyle
!!  - J_x w_{\underline{j}^{h+}}\,
!!    \frac{\partial\phi_{\underline{i}}}{\partial x_h}
!!    (\underline{x}^G_{\underline{j}^{h+}}).
!!  \end{array}
!! \f}
!! To evaluate the derivative, we have to to take into account the
!! fact that, except for the \f$h\f$-th one, all the components of the
!! two multiindexes \f$\underline{i}\f$ and \f$\underline{j}^{h+}\f$
!! have the same range. Thus, with respect to such components, there
!! are no differences with the standard DG case, and one is left with
!! \f{displaymath}{
!!    \frac{\partial\phi_{\underline{i}}}{\partial x_h}
!!    (\underline{x}^G_{\underline{j}^{h+}}) = J_h^{-1}
!!    \delta_{\underline{i}(1) \underline{j}^{h+}(1)} \ldots
!!    \phi'_{\underline{i}(h)}(x^G_{\underline{j}^{h+}(h)}) \ldots
!!    \delta_{\underline{i}(d) \underline{j}^{h+}(d)}.
!! \f}
!! This means that the volume integral is zero unless
!! \f$\underline{i}_{\backslash h}=
!! \underline{j}^{h+}_{\backslash h}\f$, in which case we have
!! \f{displaymath}{
!!   - \int_{\underline{Q}_1,K_x}
!!   \underline{q}_{h\underline{j}^{h+}}
!!   \cdot\nabla\phi_{\underline{i}} \, dx = -
!!   \frac{J_xw_{\underline{i}}}{J_hw_{\underline{i}(h)}}
!!   D^{RT}_{\underline{i}(h)\underline{j}^{h+}(h)},
!! \f}
!! with the new differentiation matrix
!! \f{displaymath}{
!!  D^{RT}_{ij}= w_j \phi_i^\prime(x^G_j), \qquad i=1,\ldots,r+1,
!!  \quad j=1,\ldots,r+2.
!! \f}
!! \note This result is formally identical to the standard DG case.
!!
!! The last integrals that must be considered are the boundary
!! integrals. Such integrals can be treated exactly in the same way as
!! for the classical DG scheme, exploiting the fact that on each side
!! \f$s\f$ both the normal component of the vector basis and the
!! scalar basis are polynomials of degree \f$r\f$.
!!
!! This leads us to the following conclusion: <em>the same
!! implementation can be used for the classical DG and the DG-RT
!! formulations, provided one introduces all the required
!! infrastructure to describe two polynomial spaces. In the case of
!! the classical DG method, such polynomial spaces coincides, while in
!! the RT case they are different</em>.
!!
!! \note We never use the fact that the vector basis is exactly one
!! degree higher that the scalar one. In fact, the two polynomial
!! degrees are arbitrary. However, it is important to use the same
!! degree for vectors and scalars in the "transversal" directions.
!!
!! \section EnergyProjection Energy conserving projection
!!
!! The electric field used in the Vlasov equation can be computed
!! using an energy conserving projection. This amounts to the
!! element-by-element computation
!! \f{displaymath}{
!!  \int_{Q_3,K_x}\underline{\tilde{E}} \cdot \underline{w}\,d x =
!!  \int_{Q_3,K_x}\nabla\Phi \cdot \underline{w}\,d x -
!!  \int_{\underline{Q}_1,K_x} \nabla\Phi
!!  \cdot\mathcal{A}\underline{w}\,d x +\int_{\underline{Q}_1,K_x}
!!  \underline{E} \cdot\mathcal{A}\underline{w}\,d x.
!! \f}
!!
!! Alternatively, one can use a simple \f$L^2\f$ projection, which
!! reduces to the identity is the extend basis coincides with the
!! scalar basis.
!!
!! Let us first consider the energy conserving projection. Taking
!! \f$\underline{w}_{k\underline{i}} =
!! \phi_{\underline{i}}\underline{e}_k\f$ we have
!! \f{displaymath}{
!!  J_xw_{\underline{i}}\tilde{E}_{k\underline{i}} = 
!!  \int_{Q_3,K_x}\frac{\partial \Phi}{\partial x_k}
!!  \phi_{\underline{i}}\,d x -
!!  \int_{\underline{Q}_1(k),K_x} \alpha_k \frac{\partial
!!  \Phi}{\partial x_k} \phi_{\underline{i}}\,d x
!!  +\int_{\underline{Q}_1(k),K_x} \alpha_k
!!  \left(\underline{E}\right)_k  \phi_{\underline{i}}\,d x.
!! \f}
!! The first term can be computed as for the matrix \f$D^{LDG}\f$, so
!! that
!! \f{displaymath}{
!!  \int_{Q_3,K_x}\frac{\partial \Phi}{\partial x_k}
!!  \phi_{\underline{i}}\,d x = \sum_{\underline{j}}
!!  \Phi_{\underline{j}}
!!  \int_{Q_3,K_x}\frac{\partial \phi_{\underline{j}}}{\partial x_k}
!!  \phi_{\underline{i}}\,d x = \frac{J_xw_{\underline{i}}}
!!  {J_kw_{\underline{i}(k)}} \sum_{h=1}^{p_k}
!!  D_{h\,\underline{i}(k)}
!!  \Phi_{\underline{i}_{\backslash k,h}},
!! \f}
!! where the summation can be seen as the matrix vector product of
!! \f$D^T\f$ and \f$\Phi\f$ along the \f$k\f$-th dimension. A similar
!! argument can be used to compute the second term, with the
!! difference that now the extended quadrature nodes are used,
!! yielding
!! \f{displaymath}{
!!  \int_{\underline{Q}_1(k),K_x}\alpha_k \frac{\partial
!!  \Phi}{\partial x_k}
!!  \phi_{\underline{i}}\,d x = \sum_{\underline{j}}
!!  \Phi_{\underline{j}}
!!  \int_{\underline{Q}_1(k),K_x}\alpha_k\frac{\partial
!!  \phi_{\underline{j}}}{\partial x_k} \phi_{\underline{i}}\,d x =
!!  \sum_{\underline{j}} \Phi_{\underline{j}}
!!  J_x \sum_{\underline{l}^{k+}} w_{\underline{l}^{k+}}\, \alpha_k(
!!  \underline{x}^G_{\underline{l}^{k+}})\,
!!  \frac{\partial \phi_{\underline{j}}}{\partial x_k}(
!!  \underline{x}^G_{\underline{l}^{k+}}) \,
!! \phi_{\underline{i}}( \underline{x}^G_{\underline{l}^{k+}}).
!! \f}
!! In the last summation, one has nonzero contribution only if
!! \f{displaymath}{
!!  \underline{i}_{\backslash k} =
!!  \underline{j}_{\backslash k} =
!!  \underline{l}^{k+}_{\backslash k}
!! \f}
!! so that
!! \f{displaymath}{
!!  \int_{\underline{Q}_1(k),K_x}\alpha_k \frac{\partial
!!  \Phi}{\partial x_k}
!!  \phi_{\underline{i}}\,d x =
!!  \frac{J_xw_{\underline{i}}} {J_kw_{\underline{i}(k)}}
!!  \sum_{s=1}^{p_k} \Phi_{\underline{i}_{\backslash k,s}}
!!  \sum_{h=1}^{p^{RT}_k} \alpha_k(
!!  \underline{x}^G_{\underline{i}_{\backslash k,h}})\,D^{RT}_{sh}
!!  P_{\underline{i}(k)\,h}
!! \f}
!! with the interpolation matrix
!! \f{displaymath}{
!!  P_{ij} = \phi_i(x^G_j), \qquad i=1,\ldots,r+1, \quad
!!  j=1,\ldots,r+2.
!! \f}
!! Now, consider that \f$\alpha_k\f$ is a function of \f$x_k\f$, so
!! that
!! \f{displaymath}{
!!  \alpha_k( \underline{x}^G_{\underline{i}_{\backslash k,h}}) =
!!  \alpha_k(x^G_h)
!! \f}
!! and, moreover, notice that the expression for \f$\alpha_k\f$ is
!! \f{displaymath}{
!!  \alpha_k(x^G_h) = \alpha^{\pm}_h = 1 \pm 2\left( x^G_h -
!!  \frac{1}{2} \right),
!! \f}
!! using the positive sign when \f$v_k>0\f$ and vice versa. Defining a
!! new differentiation matrix
!! \f{displaymath}{
!!  \mathcal{D}^{\pm}_{ij} = \sum_{h=1}^{p^{RT}_k} \alpha_h^{\pm}
!!  D^{RT}_{ih} P_{jh},
!! \f}
!! we end up with
!! \f{displaymath}{
!!  \int_{\underline{Q}_1(k),K_x}\alpha_k \frac{\partial
!!  \Phi}{\partial x_k}
!!  \phi_{\underline{i}}\,d x =
!!  \frac{J_xw_{\underline{i}}} {J_kw_{\underline{i}(k)}}
!!  \sum_{h=1}^{p_k}
!!  \mathcal{D}^{\pm}_{h\,\underline{i}(k)}
!!  \Phi_{\underline{i}_{\backslash k,h}},
!! \f}
!! which is formally identical to the expression obtained for the
!! first term. Finally, consider the third term, for which we have
!! \f{displaymath}{
!!  \int_{\underline{Q}_1(k),K_x} \alpha_k
!!  \left(\underline{E}\right)_k  \phi_{\underline{i}}\,d x =
!!  \sum_{\underline{j}^{k+}} E_{k\underline{j}^{k+}} 
!!  \int_{\underline{Q}_1(k),K_x} \alpha_k \phi_{\underline{j}^{k+}}
!!  \phi_{\underline{i}}\,d x.
!! \f}
!! After introducing the quadrature nodes of \f$\underline{Q}_1(k)\f$
!! and considering the nonzero contributions yields
!! \f{displaymath}{
!!  \int_{\underline{Q}_1(k),K_x} \alpha_k
!!  \left(\underline{E}\right)_k  \phi_{\underline{i}}\,d x =
!!  J_x \frac{w_{\underline{i}}}{w_{\underline{i}(k)}}
!!  \sum_{h=1}^{p_k^{RT}}\mathcal{P}^{\pm}_{\underline{i}(k) \, h}
!!  E_{k\underline{i}_{\backslash k,h}},
!! \f}
!! where we have defined
!! \f{displaymath}{
!!  \mathcal{P}^{\pm}_{ij} = w_j\alpha^{\pm}_jP_{ij}.
!! \f}
!! Collecting all the contributions, we obtain the local projection
!! \f{displaymath}{
!!  \tilde{E}_{k\underline{i}} = 
!!  \frac{1}{J_kw_{\underline{i}(k)}} \sum_{h=1}^{p_k}
!!  \left(D - \mathcal{D}^{\pm} \right)_{h\,\underline{i}(k)}
!!  \Phi_{\underline{i}_{\backslash k,h}}
!!  +
!!  \frac{1}{w_{\underline{i}(k)}}
!!  \sum_{h=1}^{p_k^{RT}}\mathcal{P}^{\pm}_{\underline{i}(k) \, h}
!!  E_{k\underline{i}_{\backslash k,h}}.
!! \f}
!<----------------------------------------------------------------------
module mod_poisson_pb

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 !$ use omp_lib

 !$ use mod_omp_utils, only: &
 !$   mod_omp_utils_initialized, &
 !$   detailed_timing_omp, &
 !$   omput_push_key,    &
 !$   omput_pop_key,     &
 !$   omput_start_timer, &
 !$   omput_close_timer, &
 !$   omput_write_time

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_tps_phs_grid, only: &
   mod_tps_phs_grid_initialized, &
   t_tps_phs_grid, t_tps_grid

 use mod_tps_base, only: &
   mod_tps_base_initialized, &
   t_tps_base

 use mod_f_state, only: &
   mod_f_state_initialized, &
   t_p_state, t_f_state, t_e_field

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   mpi_comm_world, mpi_comm_rank
 
 use mod_sparse, only: &
   mod_sparse_initialized, &
   ! sparse types
   t_col,       &
   t_tri,       &
   ! construction of new objects
   new_col,     &
   new_tri,     &
   ! convertions
   col2tri,     &
   tri2col,     &
   tri2col_skeleton, &
   tri2col_skeleton_part, &
   ! overloaded operators
   operator(+), &
   operator(*), &
   sum,         &
   transpose,   &
   matmul,      &
   ! error codes
   wrong_n,     &
   wrong_m,     &
   wrong_nz,    &
   wrong_dim,   &
   ! other functions
   nnz_col,     &
   nz_col,      &
   nz_col_i,    &
   get,         &
   set,         &
   diag,        &
   spdiag,      &
   ! deallocate
   clear

 use mod_octave_io_sparse, only: &
   mod_octave_io_sparse_initialized, &
   write_octave

 use mod_linsolver, only: &
   mod_linsolver_initialized, &
   c_linpb, c_itpb, c_mumpspb, c_pastixpb, c_umfpackpb, &
   gmres

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_poisson_pb_constructor, &
   mod_poisson_pb_destructor,  &
   mod_poisson_pb_initialized, &
   compute_electric_field, compute_c11, p0

 private

!-----------------------------------------------------------------------

! Module types and parameters

 !> Linear problem for the Poisson equation
 type, extends(c_itpb) :: t_poisson_it
 contains
  procedure, nopass :: pres
  procedure, nopass :: pkry
 end type t_poisson_it

 type, extends(c_mumpspb) :: t_poisson_mumps
 contains
  procedure, pass(s) :: xassign => mumps_xassign
 end type t_poisson_mumps

 type, extends(c_umfpackpb) :: t_poisson_umfpack
 contains
  procedure, pass(s) :: xassign => umf_xassign
 end type t_poisson_umfpack

 type, extends(c_pastixpb) :: t_poisson_pastix
 contains
  procedure, pass(s) :: xassign => pastix_xassign
 end type t_poisson_pastix

! Module variables

 character(len=*), parameter :: linear_solver = 'pastix'
 class(c_linpb), allocatable :: linpb
 real(wp), allocatable, target :: rhs(:)
 type(t_col), save :: mmi, dd, ddb, cc, ll, ddd, ccc, mid
 type(t_col), target, save :: mmmt
 integer, allocatable, target :: gij(:)
 real(wp), allocatable :: da(:,:,:), pa(:,:,:), dmd(:,:,:)
 !> Initial guess for the iterative solver: it is set to zero in the
 !! constructor and never changed again.
 type(t_p_state) :: p0

 ! public members
 logical, protected ::               &
   mod_poisson_pb_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_poisson_pb'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_poisson_pb_constructor(grid,base,fu,alpha_rec)
  type(t_tps_phs_grid), intent(in), target :: grid
  type(t_tps_base), intent(in) :: base !< space base
  integer, intent(in) :: fu
  logical, intent(in) :: alpha_rec

  integer :: ie, is, ivm(grid%d), ixm(grid%d), jxm(grid%d), &
    i, j, k, n, ie1, ie2, ixs, ix1, ix2, ix1x, ix2x
  integer :: ndofs_e, ndofs_p, pos
  integer, allocatable :: iim(:), jjm(:), iid(:), jjd(:), &
    iidb(:), jjdb(:), iic(:), jjc(:), iil(:), jjl(:)
  real(wp) :: xwg, c11
  real(wp), allocatable :: xxm(:), xxd(:), xxdb(:), xxc(:), xxl(:)

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
!$      ( (detailed_timing_omp.eqv..true.).and. &
!$ (mod_omp_utils_initialized.eqv..false.) ) .or. &
  (mod_state_vars_initialized.eqv..false.) .or. &
(mod_tps_phs_grid_initialized.eqv..false.) .or. &
    (mod_tps_base_initialized.eqv..false.) .or. &
     (mod_f_state_initialized.eqv..false.) .or. &
   (mod_mpi_utils_initialized.eqv..false.) .or. &
      (mod_sparse_initialized.eqv..false.) .or. &
(mod_octave_io_sparse_initialized.eqv..false.) .or. &
   (mod_linsolver_initialized.eqv..false.) ) then  
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_poisson_pb_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   allocate(rhs(base%pkd*grid%gx%ne+1))
   select case(linear_solver)
    case('gmres')
     allocate(t_poisson_it::linpb)
     select type(linpb); type is(t_poisson_it)
      linpb%abstol    = .true.
      linpb%tolerance = 1.0e-8_wp
      linpb%nmax      = 50 ! size of the Krylov space
      linpb%rmax      = 5  ! restarts
      linpb%solver    => gmres
      linpb%mpi_comm  = mpi_comm_world
      call mpi_comm_rank(linpb%mpi_comm,linpb%mpi_id,i)
     end select
    case('mumps')
     allocate(t_poisson_mumps::linpb)
     select type(linpb); type is(t_poisson_mumps)
      linpb%distributed    = .true.
      !linpb%poo           = <choose reordering>
      linpb%transposed_mat = .true.
      !linpb%gn             = mmmt%m
      !linpb%m              => mmmt
      !linpb%rhs            => rhs
      !linpb%gij            => gdofs_nat
      !linpb%gij            => gij
      !linpb%mpi_comm       = mpi_comm_world
     end select
    case('umfpack')
     allocate(t_poisson_umfpack::linpb)
     select type(linpb); type is(t_poisson_umfpack)
      linpb%print_level    = 0
      linpb%transposed_mat = .true.
      linpb%m              => mmmt
      linpb%rhs            => rhs
     end select
    case('pastix')
     allocate(t_poisson_pastix::linpb)
     select type(linpb); type is(t_poisson_pastix)
      linpb%transposed_mat = .true.
      linpb%gn             = size(rhs) ! mmmt not set yet
      linpb%m              => mmmt
      linpb%rhs            => rhs
      allocate(gij(0:linpb%gn-1)); gij = (/ (i, i=0,linpb%gn-1) /)
      linpb%gij            => gij
      linpb%mpi_comm       = mpi_comm_world
     end select
    case default
     call error(this_sub_name,this_mod_name, &
                'Unknown linear solver.')
   end select

   ! We can now define the system matrix.

   ! allocations and definitions
   ndofs_e = grid%d*base%pkdx*grid%gx%ne
   ndofs_p =        base%pkd *grid%gx%ne
   allocate(iim(ndofs_e),jjm(ndofs_e),xxm(ndofs_e))
   ! For dd it is easier to count as if the local matrices were full;
   ! then, checking whether i_\k = j_\k, some terms will be set to
   ! zero, but nevertheless for simplicity they are treated as nnz.
   i = base%pkd*grid%d*base%pkdx * grid%gx%ne ! tmp
   allocate(iid(i),jjd(i),xxd(i))
   i = 4*base%spkd * grid%gx%ns
   allocate(iic(i),jjc(i),xxc(i))
   ! notice ddb is bigger that cc, but the number of nonzero elements
   ! is the same, since n has only one nonzero component on each side
   allocate(iidb(i),jjdb(i),xxdb(i))
   i = base%pkd * grid%gx%ne
   allocate(iil(i),jjl(i),xxl(i))

   ! volume integrals
   do ie=1,grid%gx%ne
     do i=1,base%pkdx ! diagonal matrix pkdx*pkdx
       do k=1,grid%d
         xwg = grid%gx%e(ie)%vol*base%wgldx(i,k)
         pos = k + (i-1)*grid%d + (ie-1)*grid%d*base%pkdx
         iim(pos) = pos ! diagonal matrix
         jjm(pos) = pos
         xxm(pos) = 1.0_wp/xwg ! inverse mass matrix
       enddo
     enddo
     do i=1,base%pkd  ! diff. matrix: pkd*pkdx
       xwg = grid%gx%e(ie)%vol*base%wgld(i)
       ixm = base%mldx(i,:)
       do j=1,base%pkdx
         do k=1,grid%d
           jxm = base%mldxx(j,k)%mldxp
           pos = k + (j-1)*grid%d + (i-1)*grid%d*base%pkdx &
                    + (ie-1)*base%pkd*grid%d*base%pkdx
           
           ivm = jxm; ivm(k) = ixm(k) ! used to check i=j, except for k
           if(all(ivm.eq.ixm)) then
             iid(pos) = i + (ie-1)*base%pkd
             jjd(pos) = k + (j-1)*grid%d + (ie-1)*base%pkdx*grid%d
             xxd(pos) = -xwg/(grid%gx%e(ie)%bw(k)*base%wgl(ixm(k))) &
                        * base%wdphix(ixm(k),jxm(k))
           else
             iid(pos) = 1 ! avoid defining useless nonzeros
             jjd(pos) = 1
             xxd(pos) = 0.0_wp
           endif
         enddo
       enddo
       ! Now the Lagrange multiplier
       pos = i + (ie-1)*base%pkd
       iil(pos) = base%pkd * grid%gx%ne + 1
       jjl(pos) = pos
       xxl(pos) = xwg
     enddo
   enddo

   ! side integrals
   do is=1,grid%gx%ns
     n = grid%gx%s(is)%n
     ie1 = grid%gx%s(is)%ie(1)
     ie2 = grid%gx%s(is)%ie(2)
     c11 = compute_c11(is,grid%gx,base%k)
     do ixs=1,base%spkd ! side dofs
       
       ix1 = base%s_dofs(2,n,ixs) ! left and right local indexes
       ix2 = base%s_dofs(1,n,ixs)
       xwg = grid%gx%s(is)%a*base%swgld(ixs,n)

       ! now set the four terms: ix1,ix1  ix1,ix2  ix2,ix1  ix2,ix2
       pos = 1 + (ixs-1)*4 + (is-1)*4*base%spkd
       iic(pos) = ix1 + (ie1-1)*base%pkd
       jjc(pos) = ix1 + (ie1-1)*base%pkd
       xxc(pos) = c11*xwg
       pos = 2 + (ixs-1)*4 + (is-1)*4*base%spkd
       iic(pos) = ix1 + (ie1-1)*base%pkd
       jjc(pos) = ix2 + (ie2-1)*base%pkd
       xxc(pos) = -c11*xwg
       pos = 3 + (ixs-1)*4 + (is-1)*4*base%spkd
       iic(pos) = ix2 + (ie2-1)*base%pkd
       jjc(pos) = ix1 + (ie1-1)*base%pkd
       xxc(pos) = -c11*xwg
       pos = 4 + (ixs-1)*4 + (is-1)*4*base%spkd
       iic(pos) = ix2 + (ie2-1)*base%pkd
       jjc(pos) = ix2 + (ie2-1)*base%pkd
       xxc(pos) = c11*xwg

       ! same structure for ddb
       ix1x = base%s_dofsx(2,n,ixs) ! local indexes, extended basis
       ix2x = base%s_dofsx(1,n,ixs)
       pos = 1 + (ixs-1)*4 + (is-1)*4*base%spkd
       iidb(pos) = ix1 + (ie1-1)*base%pkd
       jjdb(pos) = n + (ix1x-1)*grid%d + (ie1-1)*grid%d*base%pkdx
       xxdb(pos) = 0.5_wp*xwg
       pos = 2 + (ixs-1)*4 + (is-1)*4*base%spkd
       iidb(pos) = ix1 + (ie1-1)*base%pkd
       jjdb(pos) = n + (ix2x-1)*grid%d + (ie2-1)*grid%d*base%pkdx
       xxdb(pos) = 0.5_wp*xwg
       pos = 3 + (ixs-1)*4 + (is-1)*4*base%spkd
       iidb(pos) = ix2 + (ie2-1)*base%pkd
       jjdb(pos) = n + (ix1x-1)*grid%d + (ie1-1)*grid%d*base%pkdx
       xxdb(pos) = -0.5_wp*xwg
       pos = 4 + (ixs-1)*4 + (is-1)*4*base%spkd
       iidb(pos) = ix2 + (ie2-1)*base%pkd
       jjdb(pos) = n + (ix2x-1)*grid%d + (ie2-1)*grid%d*base%pkdx
       xxdb(pos) = -0.5_wp*xwg
       
     enddo
   enddo

   ! Notice the addition of the Lagrange multiplier in dd, ddb and cc
   mmi = tri2col(new_tri(ndofs_e  ,ndofs_e  ,iim-1 ,jjm-1 ,xxm ))
   dd  = tri2col(new_tri(ndofs_p+1,ndofs_e  ,iid-1 ,jjd-1 ,xxd ))
   ddb = tri2col(new_tri(ndofs_p+1,ndofs_e  ,iidb-1,jjdb-1,xxdb))
   cc  = tri2col(new_tri(ndofs_p+1,ndofs_p+1,iic-1 ,jjc-1 ,xxc ))
   ll  = tri2col(new_tri(ndofs_p+1,ndofs_p+1,iil-1 ,jjl-1 ,xxl ))

   ! define the system matrix (exploit symmetry)
   ddd = dd + ddb
   ccc = cc + ll + transpose(ll)
   mmmt = matmul(ddd,matmul(mmi,transpose(ddd))) + ccc
   ! recover E
   mid = (-1.0_wp)*matmul(mmi,transpose(ddd))

   call write_octave(mmi ,'mmi',fu)
   call write_octave(dd  ,'dd' ,fu)
   call write_octave(ddb ,'ddb',fu)
   call write_octave(cc  ,'cc' ,fu)
   call write_octave(ll  ,'ll' ,fu)
   call write_octave(mmmt,'mmm',fu)

   deallocate(iim ,jjm ,xxm )
   deallocate(iid ,jjd ,xxd )
   deallocate(iidb,jjdb,xxdb)
   deallocate(iic ,jjc ,xxc )
   deallocate(iil ,jjl ,xxl )

   call linpb%factor('analysis')
   call linpb%factor('factorization')

   ! Compute some local matrices which will be used in the energy
   ! conserving projection
   ix1 = base%pk; ix2 = base%pkx ! used as temporaries
   allocate(da(ix1,ix1,2),pa(ix1,ix2,2),dmd(ix1,ix1,2))
   do j=1,ix1
     do i=1,ix1
       da(i,j,:) = 0.0_wp
       do k=1,ix2
         do n=1,2 ! used as temporary index
           da(i,j,n) = da(i,j,n)                                       &
             + alpha(base%xglx(k),n) * base%wdphix(i,k) * base%phix(j,k)
         enddo
       enddo
     enddo
   enddo
   if(alpha_rec) then
     dmd(:,:,1) = base%wdphi - da(:,:,1)
     dmd(:,:,2) = base%wdphi - da(:,:,2)
   else
     dmd = 0.0_wp
   endif
   do j=1,ix2
     do i=1,ix1
       do n=1,2 ! used as temporary index
         pa(i,j,n) = base%wglx(j)*alpha(base%xglx(j),n)*base%phix(i,j)
       enddo
     enddo
   enddo
   call write_octave(da ,'da' ,fu)
   call write_octave(pa ,'pa' ,fu)
   call write_octave(dmd,'dmd',fu)

   ! Currently there are no specific procedures to create and destroy
   ! t_p_state objects; one could consider adding them in mod_f_state.
   allocate(p0%pl(base%pkd*grid%gx%ne+1))
   p0%pl = 0.0_wp

   mod_poisson_pb_initialized = .true.
 contains

  pure function alpha(xi,pm) result(a)
   real(wp), intent(in) :: xi
   integer, intent(in) :: pm
   real(wp) :: a
   
    if(alpha_rec) then
      if(pm.eq.1) then
        a = 1.0_wp - 2.0_wp*(xi-0.5_wp)
      elseif(pm.eq.2) then
        a = 1.0_wp + 2.0_wp*(xi-0.5_wp)
      else
        a = -huge(1.0_wp)
      endif
    else
      a = 1.0_wp
    endif
  end function alpha

 end subroutine mod_poisson_pb_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_poisson_pb_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_poisson_pb_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate(p0%pl)

   deallocate(da,pa,dmd)

   call linpb%clean()

   deallocate(rhs,linpb)
   if(allocated(gij)) deallocate(gij)
   call clear(mmmt)
   call clear(mmi); call clear(dd); call clear(ddb)
   call clear(cc);  call clear(ll); call clear(ddd); call clear(ccc)
   call clear(mid)

   mod_poisson_pb_initialized = .false.
 end subroutine mod_poisson_pb_destructor

!-----------------------------------------------------------------------

 !> Update the electric field variables.
 !!
 !! Compute the electric field and potential as well as the particle
 !! density.
 !!
 !! \note For the solution of the linear system, the initial guess is
 !! taken from the module variable \c p0. Updating this variable must
 !! be done in some other module, typically once after each time step.
 !! This variable is never changed in this module.
 subroutine compute_electric_field(e_field,uuu)
  type(t_f_state), intent(in) :: uuu
  type(t_e_field), intent(inout) :: e_field
 
  integer :: ie, ie1, ie2, ix, iv, nxdofs, nvdofs, i, k
  integer :: k1,k2,k3 , kx1,kx2,kx3
  integer :: ixm(uuu%grid%d)
  real(wp) :: p_loc(uuu%base(1)%pkd)
  real(wp) :: vjw, xjw, cc1, cc2
  !$ integer :: dotmap_work(uuu%base(1)%d)
  character(len=*), parameter :: &
    this_sub_name = 'compute_electric_field'

   !$ if(detailed_timing_omp) then
   !$   call omput_push_key("compute_electric_fileds")
   !$   call omput_start_timer()
   !$ endif

   !$omp parallel &
   !$omp   private( nxdofs,nvdofs , ie , ie1,ie2 , iv , vjw,i,ix , &
   !$omp            xjw , p_loc , ixm , k , cc2,cc1 , k1,k2,k3 ,   &
   !$omp            kx1,kx2,kx3 , dotmap_work )                    &
   !$omp   shared( uuu , e_field , rhs,p0,linpb , mid , dmd , pa ) &
   !$omp   default(none)

   nxdofs = uuu%base(1)%pkd
   nvdofs = uuu%base(2)%pkd

   ! Diagnose the particle density: there are two equivalent
   ! implementations, the second one is parallel, while the first one
   ! would require a reduction.
   !   do ie=1,uuu%grid%ne
   !     ie1 = uuu%grid%e(ie)%e1(1)%p%o
   !     do iv=1,nvdofs
   !       vjw = uuu%grid%e(ie)%e1(2)%p%vol*uuu%base(2)%wgld(iv)
   !       e_field%r(:,ie1) = e_field%r(:,ie1) + uuu%f(:,iv,ie)*vjw
   !     enddo
   !   enddo
   !$ if(detailed_timing_omp) then
   !$omp master
   !$   call omput_push_key("diag_rho")
   !$   call omput_start_timer()
   !$omp end master
   !$ endif
   !$omp do schedule(static)
   do ie1=1,uuu%grid%gx%ne
     e_field%r(:,ie1) = 0.0_wp
     do ie2=1,uuu%grid%gv%ne
       ie = uuu%grid%i2ie( uuu%grid%gx%e(ie1)%i , uuu%grid%gv%e(ie2)%i )
       do iv=1,nvdofs
         vjw = uuu%grid%gv%e(ie2)%vol*uuu%base(2)%wgld(iv)
         e_field%r(:,ie1) = e_field%r(:,ie1) + uuu%f(:,iv,ie)*vjw
       enddo
     enddo
   enddo
   !$omp end do
   !$ if(detailed_timing_omp) then
   !$omp master
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$omp end master
   !$ endif

   ! compute the rhs: multiply by the mass matrix adding the x weights
   !$ if(detailed_timing_omp) then
   !$omp master
   !$   call omput_push_key("poisson_rhs")
   !$   call omput_start_timer()
   !$omp end master
   !$ endif
   !$omp do schedule(static)
   do ie=1,uuu%grid%gx%ne
     do ix=1,nxdofs
       xjw = uuu%grid%gx%e(ie)%vol*uuu%base(1)%wgld(ix)
       rhs(ix+(ie-1)*nxdofs) = xjw * (e_field%r(ix,ie)-1.0_wp)
     enddo
   enddo
   !$omp end do
   !$ if(detailed_timing_omp) then
   !$omp master
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$omp end master
   !$ endif

   ! the last entry of the right-hand-side corresponds to the Lagrange
   ! multiplier, so it must be set to zero
   !$omp single
   rhs(uuu%base(1)%pkd*uuu%grid%gx%ne+1) = 0.0_wp

   ! solve the linear system
   e_field%p = p0 ! initial guess, useful for iterative solvers
   call linpb%solve(e_field%p)

   ! recover the electric field
   e_field%e = reshape(matmul( mid , e_field%p%pl ) , shape(e_field%e))
   !$omp end single

   ! compute \tilde{E}
   !$ if(detailed_timing_omp) then
   !$omp master
   !$   call omput_push_key("compute_Etilde")
   !$   call omput_start_timer()
   !$omp end master
   !$ endif
   !$omp do schedule(static)
   do ie=1,uuu%grid%gx%ne
     p_loc = e_field%p%pl((ie-1)*nxdofs+1:ie*nxdofs) ! local phi dofs
     do ix=1,nxdofs
       ixm = uuu%base(1)%mldx(ix,:) ! space multiindex
       do k=1,uuu%grid%d
         cc2 = 1.0_wp/uuu%base(1)%wgl(ixm(k))
         cc1 = cc2/uuu%grid%gx%e(ie)%bw(k)
         call uuu%base(1)%dotmap (k1 ,k2 ,k3  , ixm,k &
                                 !$ , dotmap_work &
                                 )
         call uuu%base(1)%dotmapx(kx1,kx2,kx3 , ixm,k &
                                 !$ , dotmap_work &
                                 )
         do i=1,2
           e_field%et(k,i,ix,ie) =                                     &
        cc1 * dot_product(dmd(:,ixm(k),i),    p_loc(   k1:k2 :k3    )) &
      + cc2 * dot_product( pa(ixm(k),:,i),e_field%e(k,kx1:kx2:kx3,ie))
         enddo
       enddo
     enddo
   enddo
   !$omp end do
   !$ if(detailed_timing_omp) then
   !$omp master
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$omp end master
   !$ endif

   !$omp end parallel

   !$ if(detailed_timing_omp) then
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$ endif

 end subroutine compute_electric_field
 
!-----------------------------------------------------------------------

 subroutine pres(r,x)
  class(c_stv), intent(in) :: x
  class(c_stv), intent(inout) :: r

   select type(x);type is(t_p_state); select type(r);type is(t_p_state)
    r%pl = rhs - matmul(x%pl,mmmt)
   end select; end select
 end subroutine pres

!-----------------------------------------------------------------------
 
 subroutine pkry(r,x)
  class(c_stv), intent(in) :: x
  class(c_stv), intent(inout) :: r

   select type(x);type is(t_p_state);select type(r);type is(t_p_state)
    r%pl = matmul(x%pl,mmmt)
   end select; end select
 end subroutine pkry

!-----------------------------------------------------------------------

 subroutine mumps_xassign(x,s,mumps_x)
  real(wp),               intent(in)    :: mumps_x(:)
  class(t_poisson_mumps), intent(inout) :: s
  class(c_stv),           intent(inout) :: x

   select type(x); type is(t_p_state)
    x%pl = mumps_x
   end select
 end subroutine mumps_xassign

!-----------------------------------------------------------------------

 subroutine umf_xassign(x,s,umfpack_x)
  real(wp),                 intent(in)    :: umfpack_x(:)
  class(t_poisson_umfpack), intent(inout) :: s
  class(c_stv),             intent(inout) :: x

   select type(x); type is(t_p_state)
    x%pl = umfpack_x
   end select
 end subroutine umf_xassign

!-----------------------------------------------------------------------

 subroutine pastix_xassign(x,s,pastix_x)
  real(wp),                 intent(in)    :: pastix_x(:)
  class(t_poisson_pastix), intent(inout) :: s
  class(c_stv),             intent(inout) :: x

   select type(x); type is(t_p_state)
    x%pl = pastix_x
   end select
 end subroutine pastix_xassign

!-----------------------------------------------------------------------

 !> Compute \f$c_{11}=\frac{c}{h}\f$.
 !!
 !! The main purpose of this function is ensuring that \f$c_{11}\f$ is
 !! always computed in the same way.
 !!
 !! In the present version, the element size \f$h\f$ is measured in
 !! the direction orthogonal to the selected size, and we set
 !! \f$c=(k+1)^2\f$, following [Ayuso, Hajian].
 !!
 !! \note Numerical evidence indicates that the choice \f$c_{11}=0\f$
 !! is possible when using a Raviart-Thomas basis, which with the
 !! usual DG basis it leads to a singular problem (however, in 1D, the
 !! problem seems to be well posed in the case of odd numbers of
 !! elements).
 pure function compute_c11(is,grid,k) result(c11)
  integer, intent(in) :: is, k
  type(t_tps_grid), intent(in) :: grid
  real(wp) :: c11

  real(wp), parameter :: c = 1.0_wp
  integer :: n
  character(len=*), parameter :: &
    this_fun_name = 'compute_c11'

   n = grid%s(is)%n
   c11 = 0.5_wp*(c/grid%s(is)%e(1)%p%bw(n) + c/grid%s(is)%e(2)%p%bw(n))
   !c11 = 1.0_wp ! independent from h
   !c11 = 0.0_wp ! don't use c11: possible with the RT space
   ! Correction for the polynomial order
   c11 = real( (k+1)**2 ,wp) * c11
 
 end function compute_c11
 
!-----------------------------------------------------------------------

end module mod_poisson_pb

